library(tidyverse)
library(broom)
library(haven)
library(stargazer)
library(kableExtra)

data2 <- read_dta("Study_2.dta")

# Create treatment variable
data2 <- data2 %>%
  mutate(treatment = case_when(
    trustexperimentsplit == 1 ~ 0,
    trustexperimentsplit == 2 ~ 1
  )) %>%
  mutate(treatment = factor(treatment, labels = c("Control", "Treatment")))

# Create trust variables
data2 <- data2 %>%
  mutate(
    trustMP = coalesce(trustExperimentGroup1MP, trustExperimentGroup2MP),
    trustMSP = coalesce(trustExperimentGroup1MSP, trustExperimentGroup2MSP),
    trustScotMin = coalesce(trustExperimentGroup1ScotMin, trustExperimentGroup2ScotMin),
    trustUKMin = coalesce(trustExperimentGroup1UKMin, trustExperimentGroup2UKMin),
    trustCivWest = coalesce(trustExperimentGroup1CivWest, trustExperimentGroup2CivWest),
    trustCivScot = coalesce(trustExperimentGroup1CivScot, trustExperimentGroup2CivScot)
  ) %>%
  mutate(across(c(trustMP, trustMSP, trustScotMin, trustUKMin, trustCivWest, trustCivScot), 
                ~if_else(. == 8, NA_real_, . - 1)))

# Prepare data for plotting
plot_data2 <- data2 %>%
  pivot_longer(cols = c(trustMP, trustMSP, trustUKMin, trustScotMin, trustCivWest, trustCivScot),
               names_to = "outcome", values_to = "trust") %>%
  mutate(outcome = factor(outcome, 
                          levels = c("trustMP", "trustMSP", "trustUKMin", "trustScotMin", "trustCivWest", "trustCivScot"),
                          labels = c("MPs", "MSPs", "UK Ministers", "Scottish Ministers", "UK Civil Servants", "Scottish Civil Servants"))) %>%
  mutate(outcome = fct_rev(outcome))

# Calculate means and confidence intervals
means_data2 <- plot_data2 %>%
  group_by(outcome, treatment) %>%
  summarise(
    mean_trust = mean(trust, na.rm = TRUE),
    se_trust = sd(trust, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean_trust - 1.96 * se_trust,
    ci_upper = mean_trust + 1.96 * se_trust,
    .groups = "drop"
  )

# Function to perform t-test and return p-value
t_test_p_value <- function(outcome_var, data) {
  t_test_result <- t.test(outcome_var ~ treatment, data = data)
  return(t_test_result$p.value)
}

# Calculate ATE and p-values
ate_data2 <- plot_data2 %>%
  group_by(outcome) %>%
  summarise(
    p_value = t_test_p_value(trust, cur_data()),
    .groups = "drop"
  ) %>%
  left_join(
    means_data2 %>%
      pivot_wider(id_cols = outcome, names_from = treatment, values_from = mean_trust),
    by = "outcome"
  ) %>%
  mutate(
    ate = Treatment - Control,
    stars = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01 ~ "**",
      p_value < 0.05 ~ "*",
      TRUE ~ ""
    ),
    ate_label = sprintf("ATE = %.2f%s", ate, stars)
  )

# Create main plot
plot2 <- ggplot() +
  geom_point(data = plot_data2, aes(x = outcome, y = trust, color = treatment), 
             position = position_jitterdodge(jitter.width = 0.3, jitter.height = 0.35, dodge.width = 0.8),
             alpha = 0.075, size = 1.5, shape = 16) +
  geom_pointrange(data = means_data2, 
                  aes(x = outcome, y = mean_trust, ymin = ci_lower, ymax = ci_upper, group = treatment, color = treatment),
                  position = position_dodge(width = 0.8),
                  size = 0.8) +
  geom_text(data = ate_data2, aes(x = outcome, y = 7.2, label = ate_label),
            vjust = -0.5, size = 3) +
  scale_color_manual(values = c("#11CE9E", "#CE1141FF"), name = "Group") +
  labs(y = "Trust score",
       caption = "* p<0.05, ** p<0.01, *** p<0.001") +
  theme_minimal() +
  theme(
    axis.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    axis.text.x = element_text(angle = 0, hjust = 1),
    legend.position = "bottom"
  ) +
  coord_flip(ylim = c(0, 7.5))

ggsave("Figure_2b.png", plot2, width = 12, height = 6)

# Create bar chart for control group distributions
control_data2 <- plot_data2 %>% 
  filter(treatment == "Control")

control_dist_plot2 <- ggplot(control_data2, aes(x = trust, y = ..prop.., group = outcome)) +
  geom_bar(fill = "#11CE9E", alpha = 0.5, color = "#11CE9E", size = 0.5) +
  facet_wrap(~outcome, ncol = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = "Trust Score (Higher = More Trust)", y = "Percentage") +
  theme_minimal() +
  theme(legend.position = "none")

ggsave("Figure_7.png", 
       control_dist_plot2, width = 10, height = 8)

# Function to perform t-test and return tidy results with group means
run_t_test <- function(data, outcome_var) {
  t_test_result <- t.test(data[[outcome_var]] ~ data$treatment)
  control_mean <- mean(data[[outcome_var]][data$treatment == "Control"], na.rm = TRUE)
  treatment_mean <- mean(data[[outcome_var]][data$treatment == "Treatment"], na.rm = TRUE)
  n <- sum(!is.na(data[[outcome_var]]))
  
  tibble(
    outcome = outcome_var,
    mean_diff = t_test_result$estimate[2] - t_test_result$estimate[1],
    control_mean = control_mean,
    treatment_mean = treatment_mean,
    p_value = t_test_result$p.value,
    n = n
  )
}

# List of outcome variables
outcome_vars <- c("trustMP", "trustMSP", "trustUKMin", "trustScotMin", "trustCivWest", "trustCivScot")

# Run t-tests for all outcome variables
results2 <- map_df(outcome_vars, ~run_t_test(data2, .x))

# Create a formatted table
table2 <- results2 %>%
  mutate(
    outcome = case_when(
      outcome == "trustMP" ~ "MPs",
      outcome == "trustMSP" ~ "MSPs",
      outcome == "trustUKMin" ~ "UK Ministers",
      outcome == "trustScotMin" ~ "Scottish Ministers",
      outcome == "trustCivWest" ~ "UK Civil Servants",
      outcome == "trustCivScot" ~ "Scottish Civil Servants"
    ),
    significance = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01 ~ "**",
      p_value < 0.05 ~ "*",
      p_value < 0.1 ~ ".",
      TRUE ~ ""
    )
  ) %>%
  select(outcome, mean_diff, control_mean, treatment_mean, p_value, significance, n) %>%
  mutate(across(where(is.numeric), ~round(., 3)))

# Create and save LaTeX table
latex_table2 <- kable(table2, format = "latex", booktabs = TRUE,
                      col.names = c("Outcome", "Mean Difference", "Control Mean", "Treatment Mean", "p-value", "", "N"),
                      caption = "Treatment Effects on Trust in Political Actors - Study 2") %>%
  kable_styling(latex_options = c("striped", "hold_position"))

# Save the table to a file
writeLines(latex_table2, "Table_3.tex")

# Print the table in the console
print(latex_table2)